package scu.maqiang.numeric;

import static java.lang.Math.sqrt;

/**
 * 复稀疏矩阵迭代求解器
 * 
 * @author 马强
 *
 */
public class IterSCSolver extends Solver {
	public static void main(String[] args) {
		/*
		 * SCMatrix mat = new SCMatrix(2); mat.setElement(0, 0, new Complex(4.0, 0.0));
		 * mat.setElement(1, 1, new Complex(4.0, 0.0)); mat.setElement(0, 1, new
		 * Complex(1.0, 2.0)); mat.setElement(1, 0, new Complex(1.0, -2.0)); mat.sort();
		 * System.out.println(mat.display()); Complex[] b = {new Complex(2, 1), new
		 * Complex(1, 2)}; Complex[] x = {new Complex(0, 0), new Complex(0, 0)};
		 * IterSCSolver solver = new IterSCSolver(mat); solver.PCGSSOR(b, x, 1.5, 1);
		 * System.out.println(Complex.toString(x));
		 */
		/*
		SCMatrix mat = new SCMatrix(2);
		mat.setElement(0, 0, new Complex(4.0, 0.0));
		mat.setElement(1, 1, new Complex(4.0, 0.0));
		mat.setElement(0, 1, new Complex(1.0, 2.0));
		mat.setElement(1, 0, new Complex(1.0, -1.0));
		mat.sort();
		System.out.println(mat.display());
		Complex[] b = { new Complex(2, 1), new Complex(1, 3) };
		Complex[] x = { new Complex(0, 0), new Complex(0, 0) };
		IterSCSolver solver = new IterSCSolver(mat);
		solver.PCGSSSOR(b, x, 1.5, 1);
		System.out.println(Complex.toString(x));
		*/
		
		SCMatrix mat = new SCMatrix(10);
		for(int i = 0; i < 10; i++) {
			mat.setElement(i, i, new Complex(4.0, 0.0));
			if(i > 0) {
				mat.setElement(i, i - 1, new Complex(1, -1));
			}
			if(i + 1 < 10) {
				mat.setElement(i, i + 1, new Complex(1, 1));
			}
		}
		/*
		Complex[] b = new Complex[10];
		b[0] = new Complex(5, 1);
		for(int i = 1; i <= 8; i++) {
			b[i] = new Complex(6, 0);
		}
		b[9] = new Complex(5, -1);
		*/
		Complex[] b = new Complex[10];
		b[0] = new Complex(-1, 5);
		for(int i = 1; i <= 8; i++) {
			b[i] = new Complex(0, 6);
		}
		b[9] = new Complex(1, 5);
		mat.sort();
		System.out.println(mat.display());
		IterSCSolver solver = new IterSCSolver(mat);
		Complex[] x = Complex.zeros(10);
		System.out.println(Complex.toString(b));
		solver.PCGSSSOR(b, x, 1.5, 1);
		//solver.CG(b, x, 1);
		System.out.println(Complex.toString(x));
		//System.out.println(mat.display());
		//System.out.println(Complex.toString((mat.mul(x))));
	}

	public IterSCSolver() {
		super();
	}
	/**
	 * 构造函数, 传入要求解的复稀疏矩阵
	 * @param mat SCMatrix对象
	 */
	public IterSCSolver(SCMatrix mat) {
		int n = mat.getN();
		diag = Complex.zeros(n);
		matU = new CSRCMatrix(n);
		matL = new CSRCMatrix(n);
		CSRCMatrix.SCMatrixToCSRCMatrix(mat, diag, matU, matL);
		ABMatrix = new SRMatrix(n * 2);
		for(int i = 0; i < n; i++) {
			for(SCEntry ele: mat.mat[i]) {
				ABMatrix.setElementFast(i, ele.col, ele.val.real);
				ABMatrix.setElementFast(i, ele.col + n, -ele.val.image);
				ABMatrix.setElementFast(i + n, ele.col, ele.val.image);
				ABMatrix.setElementFast(i + n, ele.col + n, ele.val.real);
			}
		}
		//ABMatrix.save("ABMatrix.txt");
	}

	
	/**
	 * 
	 * @param mat
	 * @param isDelete
	 */
	public IterSCSolver(SCMatrix mat, boolean isDelete) {
		this(mat);
		if (isDelete) {
			mat.clear();
		}
	}

	/**
	 * ���캯�������������ص�ϡ�������������Լ�����������
	 * 
	 * @param mat          SCMatrix����
	 * @param eps          �������
	 * @param maxIteration ����������
	 */
	public IterSCSolver(SCMatrix mat, double eps, int maxIteration) {
		this(mat);
		this.eps = eps;
		this.maxIteration = maxIteration;
	}

	public void setMatrix(SCMatrix mat) {
		int n = mat.getN();
		if (diag != null) {
			for(int i = 0, ni = diag.length; i < ni; i++) {
				diag[i] = null;
			}
		}
		diag = Complex.zeros(n);
		if(matU != null) {
			matU.clear();
		}
		matU = new CSRCMatrix(n);
		if (matL != null) {
			matL.clear();
		}
		matL = new CSRCMatrix(n);
		CSRCMatrix.SCMatrixToCSRCMatrix(mat, diag, matU, matL);
	}

	public int Jacobi(SCMatrix A, Complex[] b, Complex[] x0) {
		int n = A.getN();
		Complex[] x1 = Complex.zeros(n);
		double err = eps + 1.0;
		Complex a00 = new Complex(0.0, 0.0);
		int i;
		for (i = 0; i < maxIteration; i++) {
			if(err < Constants.Er) {
				break;
			}
			for(int j = 0; j < n; j++) {
				Complex temp = new Complex(0.0, 0.0);
				for (SCEntry ele: A.mat[j]) {
					if(ele.col != j) {
						temp.real += ele.val.real * x0[ele.col].real - ele.val.image * x0[ele.col].image;
						temp.image += ele.val.real * x0[ele.col].image + ele.val.image *x0[ele.col].real;
					}else {
						a00 = ele.val;
					}
				}
				Complex diff = new Complex(b[j].real - temp.real, b[j].image - temp.image);
				x1[j] = diff.div(a00);
			}
			err = Complex.norm(Complex.sub(x0, x1));
			//System.out.println(x0.length + "\t" +x0.length + "\t" + n + "\n");
			for(int j = 0; j < n; j++) {
				x1[j].real = x0[j].real;
				x1[j].image = x0[j].image;
			}
		}
		if(i == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time! ");
			return -1;
		}
		System.out.println("Iteration Time " + i);
		return i;
	}
	
	/**
	 * �����ݶȷ����Գ��������Է�����(�Ѵ�����)
	 * 
	 * @param b         �Ҷ���
	 * @param x         ������ʼֵ, ͬʱ���������
	 * @param verbosity ���������
	 * @return ��������, ��������������������������-1
	 */
	public void CG(Complex[] b, Complex[] x, int verbosity) {
		long begin = System.nanoTime();
		int n = matU.getM();
		Complex[] Ap = matU.mul(x);
		Complex[] r = new Complex[n];
		Complex[] p = new Complex[n];
		for (int i = 0; i < n; i++) {
			r[i] = new Complex(b[i].real - Ap[i].real, b[i].image - Ap[i].image);
			p[i] = new Complex(r[i]);
		}
		double r2 = Complex.dotT(r, r).real;
		double temp1, alpha, beta;
		int i, k;
		for (i = 0; i < maxIteration; i++) {
			if (sqrt(r2) < eps) {
				break;
			}
			matU.mul(p, Ap);
			temp1 = Complex.dotT(Ap, p).real;
			alpha = r2 / temp1;
			temp1 = 0;
			for (k = 0; k < n; k++) {
				x[k].real += alpha * p[k].real;
				x[k].image += alpha * p[k].image;
				r[k].real -= alpha * Ap[k].real;
				r[k].image -= alpha * Ap[k].image;
				temp1 += r[k].mulT(r[k]).real;
			}
			beta = temp1 / r2;
			r2 = temp1;
			for (k = 0; k < n; k++) {
				p[k].real = r[k].real + beta * p[k].real;
				p[k].image = r[k].image + beta * p[k].image;
			}
			if (verbosity >= 2) {
				System.out.println();
				System.out.println("Iteration: " + (i + 1));
				System.out.printf("min: %20.10f\t max: %20.10f\n", Complex.min(x), Complex.max(x));
				System.out.printf("Error: %20.10f\n", sqrt(r2));
				System.out.println();
			}
		}
		if (verbosity >= 1) {
			System.out.println();
			System.out.println("Total iteration: " + i);
			System.out.printf("min: %20.10f\t max: %20.10f\n", Complex.min(x), Complex.max(x));
			long end = System.nanoTime();
			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
			System.out.println();
		}
	}

	/**
	 * �����ݶȷ��������(CGNE)�����һ�����Է�����(�Ѵ�����)
	 * 
	 * @param b         �Ҷ���
	 * @param x0        ������ʼֵ�������ص�����
	 * @param verbosity ������Ʋ���
	 * @return ��������, ��������������������������-1
	 */
	public int CGNE(Complex[] b, Complex[] x0, int verbosity) {
		int n = b.length;
		long begin = System.nanoTime();
		Complex[] Ap0 = matU.mul(x0);
		Complex[] r0 = new Complex[n];
		for (int i = 0; i < n; i++) {
			r0[i] = new Complex(b[i].real - Ap0[i].real, b[i].image - Ap0[i].image);
		}
		Complex[] p0 = matU.mulT(r0);
		double temp1, temp2, alpha, beta, r2;
		r2 = Complex.dotT(r0, r0).real;
		int i;
		for (i = 0; i < maxIteration; i++) {
			if (sqrt(r2)/Math.sqrt(Complex.dotT(b, b).real) < eps) {
				break;
			}
			temp1 = 0.0;
			temp2 = 0.0;
			matU.mul(p0, Ap0);
			temp1 = Complex.dotT(r0, r0).real;
			temp2 = Complex.dotT(p0, p0).real;
			alpha = temp1 / temp2;
			temp2 = temp1;
			temp1 = 0.0;
			for (int j = 0; j < n; j++) {
				x0[j].real += alpha * p0[j].real;
				x0[j].image += alpha * p0[j].image;
				r0[j].real -= alpha * Ap0[j].real;
				r0[j].image -= alpha * Ap0[j].image;
				temp1 += r0[j].mulT(r0[j]).real;
			}
			beta = temp1 / temp2;
			matU.mulT(r0, Ap0);
			for (int j = 0; j < n; j++) {
				p0[j].real = Ap0[j].real + beta * p0[j].real;
				p0[j].image = Ap0[j].image + beta * p0[j].image;
			}
			r2 = Complex.dotT(r0, r0).real;
			if (verbosity >= 2) {
				System.out.println();
				System.out.println("Iteration: " + (i + 1));
				System.out.printf("min: %20.10e\t max: %20.10e\n", Complex.min(x0), Complex.max(x0));
				System.out.printf("Error: %20.10e\n", r2);
				System.out.println();
			}
		}
		if (verbosity >= 1) {
			System.out.println();
			System.out.println("Total iteration: " + i);
			System.out.printf("min: %20.10e\t max: %20.10e\n", Complex.min(x0), Complex.max(x0));
			long end = System.nanoTime();
			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
			System.out.println();
		}
		if (i == maxIteration) {
			throw new IllegalArgumentException("Iteration Has Reached Maximum Time");
		}
		return i;
	}

	/**
	 * 
	 * @param  b 
	 * @param x0 
	 * @param omega
	 * @param verbosity
	 * @return 
	 */
	public int PCGSSOR(Complex[] b, Complex[] x0, double omega, int verbosity) {
		int n = matU.getN();
		long begin = System.nanoTime();
		Complex[] r0 = new Complex[n];
		Complex[] z0 = matU.mul(x0);
		Complex[] z1 = new Complex[n];
		Complex[] p0 = new Complex[n];
		for (int i = 0; i < n; i++) {
			r0[i] = new Complex(b[i].real - z0[i].real, b[i].image - z0[i].image);
			z0[i].real = 0.0;
			z0[i].image = 0.0;
			z1[i] = new Complex(0.0, 0.0);
			p0[i] = new Complex(0.0, 0.0);
		}

		CSRCMatrix.SSORTreat(diag, matU, matL, r0, z0, omega);
		for (int i = 0; i < n; i++) {
			z1[i].real = z0[i].real;
			z1[i].image = z0[i].image;
			p0[i].real = z0[i].real;
			p0[i].image = z0[i].image;
		}
		int i, j;
		double temp1, temp2, alpha, beta;
		for (i = 0; i < maxIteration; i++) {
			if (Math.sqrt(Complex.dotT(r0, r0).real) < eps) {
				break;
			}
			temp1 = temp2 = 0.0;
			matU.mul(p0, z1);
			temp1 = Complex.dotT(z0, r0).real;
			temp2 = Complex.dotT(z1, p0).real;
			alpha = temp1 / temp2;
			for (j = 0; j < n; j++) {
				x0[j].real += alpha * p0[j].real;
				x0[j].image += alpha * p0[j].image;
				r0[j].real -= alpha * z1[j].real;
				r0[j].image -= alpha * z1[j].image;
			}

			// System.out.println("r0 = " + MVO.toString(r0));
			// SSORTreat(A, diag, r0, z1, omega);
			CSRCMatrix.SSORTreat(diag, matU, matL, r0, z1, omega);

			temp2 = temp1;
			temp1 = Complex.dotT(z1, r0).real;
			beta = temp1 / temp2;
			for (j = 0; j < n; j++) {
				p0[j].real = z1[j].real + beta * p0[j].real;
				p0[j].image = z1[j].image + beta * p0[j].image;
				z0[j].real = z1[j].real;
				z0[j].image = z1[j].image;
			}
			if (verbosity >= 2) {
				System.out.println();
				System.out.println("Iteration: " + (i + 1));
				System.out.printf("min: %20.10f\t max: %20.10f\n", Complex.min(x0), Complex.max(x0));
				System.out.printf("Error: %20.10f\n", Math.sqrt(Complex.dotT(r0, r0).real));
				System.out.println();
			}
			// System.arraycopy(r1, 0, r0, 0, n);
			// System.arraycopy(z1, 0, z0, 0, n);
		}

		if (verbosity >= 1) {
			System.out.println();
			System.out.println("Total iteration: " + i);
			System.out.printf("min: %20.10f\t max: %20.10f\n", Complex.min(x0), Complex.max(x0));
			long end = System.nanoTime();
			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
			System.out.println();
		}
		if (i == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time! ");
			return -1;
		}
		r0 = null;
		z0 = null;
		z1 = null;
		p0 = null;
		return i;
	}

	/**
	 * �Գ���γ��ɳ�Ԥ�������ݶ�ƽ����(CGS)���һ�����Է�����
	 *
	 * @param b  �Ҷ���
	 * @param x0 ������ʼֵ�������ص�����
	 * @param omega
	 * @param verbosity
	 * @return ��������, ��������������������������-1
	 */
	public int PCGSSSOR(Complex[] b, Complex[] x0, double omega, int verbosity) {
		int n = matU.getN();
		long begin = System.nanoTime();
		Complex[] r0 = new Complex[n];
		Complex[] z0 = matU.mul(x0);
		Complex[] u0 = new Complex[n];
		Complex[] p0 = new Complex[n];
		Complex tt = new Complex(0, 0);
		for (int i = 0; i < n; i++) {
			r0[i] = new Complex(b[i].real - z0[i].real, b[i].image - z0[i].image);
			z0[i].real = 0.0;
			z0[i].image = 0.0;
			u0[i] = new Complex(0.0, 0.0);
			p0[i] = new Complex(0.0, 0.0);
		}

		CSRCMatrix.SSORTreat(diag, matU, matL, r0, z0, omega);

		for (int i = 0; i < n; i++) {
			u0[i].real = z0[i].real;
			u0[i].image = z0[i].image;
			p0[i].real = z0[i].real;
			p0[i].image = z0[i].image;
		}
		Complex[] Ap0 = matU.mul(p0);

		Complex[] q0 = new Complex[n];
		Complex[] v0 = new Complex[n];
		Complex[] rr0 = new Complex[n];
		for (int i = 0; i < n; i++) {
			q0[i] = new Complex(0.0, 0.0);
			v0[i] = new Complex(0.0, 0.0);
			rr0[i] = new Complex(z0[i].real + 1.0, z0[i].image + 1.0);
		}

		CSRCMatrix.SSORTreat(diag, matU, matL, Ap0, v0, omega);
		int k;
		double r0Norm, bNorm;
		double temp1, temp2, alpha, beta;
		for (k = 0; k < maxIteration; k++) {
			Complex.dotT(r0, r0, tt);
			r0Norm = Math.sqrt(tt.real);
			Complex.dotT(b, b, tt);
			bNorm = Math.sqrt(tt.real);
			if (r0Norm / bNorm < eps) {
				break;
			}
			Complex.dotT(z0, rr0, tt);
			temp1 = tt.real;
			Complex.dotT(v0, rr0, tt);
			temp2 = tt.real;
			alpha = temp1 / temp2;

			for (int i = 0; i < n; i++) {
				q0[i].real = u0[i].real - alpha * v0[i].real;
				q0[i].image = u0[i].image - alpha * v0[i].image;
				v0[i].real = u0[i].real + q0[i].real;
				v0[i].image = u0[i].image + q0[i].image;
			}

			matU.mul(v0, u0);
			for (int i = 0; i < n; i++) {
				x0[i].real += alpha * v0[i].real;
				x0[i].image += alpha * v0[i].image;
				r0[i].real -= alpha * u0[i].real;
				r0[i].image -= alpha * u0[i].image;
			}
			CSRCMatrix.SSORTreat(diag, matU, matL, r0, z0, omega);

			temp2 = temp1;
			Complex.dotT(z0, rr0, tt);
			temp1 = tt.real;
			beta = temp1 / temp2;
			for (int i = 0; i < n; i++) {
				u0[i].real = z0[i].real + beta * q0[i].real;
				u0[i].image = z0[i].image + beta * q0[i].image;
				p0[i].real = u0[i].real + beta * (q0[i].real + beta * p0[i].real);
				p0[i].image = u0[i].image + beta * (q0[i].image + beta * p0[i].image);
			}

			matU.mul(p0, Ap0);
			CSRCMatrix.SSORTreat(diag, matU, matL, Ap0, v0, omega);
			if (verbosity >= 2) {
				System.out.println();
				System.out.println("Iteration: " + (k + 1));
				System.out.printf("min: %20.10e\t max: %20.10e\n", Complex.min(x0), Complex.max(x0));
				System.out.printf("Error: %20.10e\n", r0Norm / bNorm);
				System.out.println();
			}
		}

		if (verbosity >= 1) {
			System.out.println();
			System.out.println("Total iteration: " + k);
			System.out.printf("min: %20.10e\t max: %20.10e\n", Complex.min(x0), Complex.max(x0));
			long end = System.nanoTime();
			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
			System.out.println();
		}

		if (k == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time! ");
			return -1;
		}
		for(int i = 0; i < n; i++) {
			r0[i] = null;
			z0[i] = null;
			u0[i] = null;
			p0[i] = null;
			Ap0[i] = null;
			q0[i] = null;
			v0[i] = null;
			rr0[i] = null;
		}
		return k;
	}

//	/**
//	 * 
//	 * @param b
//	 * @param x
//	 * @param m
//	 */
//	public void GMRES(double[] b, double[] x, int m) {
//		int n = matU.getM();
//		ParamCheck.checkEqual(n, b.length);
//		ParamCheck.checkEqual(n, x.length);
//		double[][] H = new double[m + 1][m];
//		double[][] Q = new double[n][m + 1];
//		double[] x0 = x.clone();
//		double[] f = new double[m + 1];
//		boolean flag = true;
//		int iter = 0;
//		while (flag && iter < maxIteration) {
//			double[] r = MVO.add(b, -1.0, matU.mul(x));
//			f[0] = MVO.L2Norm(r);
//			for (int i = 0; i < n; i++) {
//				Q[i][0] = r[i] / f[0];
//			}
//			for (int k = 0; k < m; k++) {
//				double[] colK = MVO.col(Q, k);
//				double[] y = matU.mul(colK);
//				for (int j = 0; j <= k; j++) {
//					double[] colJ = MVO.col(Q, j);
//					H[j][k] = MVO.dot_product(colJ, y);
//					for (int t = 0; t < n; t++) {
//						y[t] -= H[j][k] * colJ[t];
//					}
//				}
//				H[k + 1][k] = MVO.L2Norm(y);
//				for (int t = 0; t < n; t++) {
//					Q[t][k + 1] = y[t] / H[k + 1][k];
//				}
//				double[] c = new double[k + 1];
//				double error = IterDSolver.MLS_QR_Givens(MVO.subMatrix(H, 0, k + 1, 0, k), MVO.subVector(f, 0, k + 1), c);
//				MVO.add(x0, 1.0, MVO.matmul(MVO.subMatrix(Q, 0, n - 1, 0, k), c), x);
//				System.out.println("Error = " + error + "    Iter = " + iter);
//				if (Math.abs(error) < Constants.Er) {
//					flag = false;
//					break;
//				}
//				iter++;
//			}
//			for(int i = 0; i < n; i++) {
//				x0[i] = x[i];
//			}
//		}
//	}
//
//	/**
//	 * 
//	 * @param b
//	 * @param x
//	 * @param omega
//	 * @param m
//	 * @param verbosity
//	 */
//	public void GMRES_SSOR(double[] b, double[] x, double omega, int m, int verbosity) {
//		int n = matU.getM();
//		ParamCheck.checkEqual(n, b.length);
//		ParamCheck.checkEqual(n, x.length);
//		long begin = System.nanoTime();
//		double[][] H = new double[m + 1][m];
//		double[][] Q = new double[n][m + 1];
//		double[] x0 = x.clone();
//		double[] f = new double[m + 1];
//		double[] z0 = new double[n];
//		double[] r = new double[n];
//		boolean flag = true;
//		int iter = 0;
//		while (flag && iter < maxIteration) {
//			MVO.add(b, -1.0, matU.mul(x), r);
//			CSCMatrix.SSORTreat(diag, matU, matL, r, z0, omega);
//			r = z0.clone();
//			f[0] = MVO.L2Norm(r);
//			for (int i = 0; i < n; i++) {
//				Q[i][0] = r[i] / f[0];
//			}
//			for (int k = 0; k < m; k++) {
//				double[] colK = MVO.col(Q, k);
//				CSCMatrix.SSORTreat(diag, matU, matL, matU.mul(colK), z0, omega);
//				double[] y = z0.clone();
//				for (int j = 0; j <= k; j++) {
//					double[] colJ = MVO.col(Q, j);
//					H[j][k] = MVO.dot_product(colJ, y);
//					for (int t = 0; t < n; t++) {
//						y[t] -= H[j][k] * colJ[t];
//					}
//				}
//
//				H[k + 1][k] = MVO.L2Norm(y);
//				for (int t = 0; t < n; t++) {
//					Q[t][k + 1] = y[t] / H[k + 1][k];
//				}
//				//System.out.println(MVO.toString(H));
//				double[] c = new double[k + 1];
//				double error = IterDSolver.MLS_QR_Givens(MVO.subMatrix(H, 0, k + 1, 0, k), MVO.subVector(f, 0, k + 1), c);
//				//double error = IterDSolver.MLS_QR_Givens2(H, 0, k + 1, 0, k, f, c);
//				//System.out.println(Arrays.toString(c));
//				MVO.add(x0, 1.0, MVO.matmul(MVO.subMatrix(Q, 0, n - 1, 0, k), c), x);
//				if (verbosity >= 2) {
//					System.out.println("Error = " + error + "    iter = " + iter);
//					System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(x), MVO.max(x));
////					if (k == 3) {
////						System.exit(0);
////					}
//				}
//				if (Math.abs(error) < eps) {
//					flag = false;
//					break;
//				}
//				iter++;
//			}
//			for(int i = 0; i < n; i++) {
//				x0[i] = x[i];
//			}
//		}
//		
//		if (verbosity >= 1) {
//			System.out.println();
//			System.out.println("Total iteration: " + iter);
//			System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(x), MVO.max(x));
//			long end = System.nanoTime();
//			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
//			System.out.println();
//		}
//	}
//	
//	/**
//	 * �����ݶȷ����̲���(CGNR)�����һ�����Է�����(�Ѵ�����)
//	 * @param A һ�������ϡ�����
//	 * @param b �Ҷ���
//	 * @param x0 ������ʼֵ�������ص�����
//	 * @return ��������, ��������������������������-1
//	 */
//	public int CGNR(double[] b, double[] x0, int verbosity) {
//		int n = b.length;
//		double[] r0 = MVO.add(b, -1.0, matU.mul(x0));
//		double[] p0 = matU.mulT(r0);
//		double[] Ap0 = new double[n];
//		int i;
//		long begin = System.nanoTime();
//		for (i = 0; i < maxIteration; i++) {
//			if (MVO.L2Norm(r0) < eps) {
//				break;
//			}
//			double temp1 = 0.0;
//			double temp2 = 0.0;
//			Arrays.fill(Ap0, 0.0);
//			matU.mulT(r0, Ap0);
//			temp1 = MVO.dot_product(Ap0, Ap0);
//			matU.mul(p0, Ap0);
//			temp2 = MVO.dot_product(Ap0, Ap0);
//			double alpha = temp1 / temp2;
//			for (int j = 0; j < n; j++) {
//				x0[j] += alpha * p0[j];
//				r0[j] -= alpha * Ap0[j];
//			}
//			
//			temp2 = temp1;
//			temp1 = 0.0;
//			//Arrays.fill(Ap0, 0.0);
//			matU.mulT(r0, Ap0);
//			temp1 = MVO.dot_product(Ap0, Ap0);
//			double beta = temp1 / temp2;
//			for (int j = 0; j < n; j++) {
//				p0[j] = Ap0[j] + beta * p0[j];
//			}
//
//			if (verbosity >= 2) {
//				System.out.println("Error = " + r0 + "    iter = " + i);
//				System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(x0), MVO.max(x0));
////				if (k == 3) {
////					System.exit(0);
////				}
//			}
////			System.out.println("Iteration " + (i + 1) + ": min\t"
////					+ MVO.min(x0) + "\tmax\t" + MVO.max(x0));
//		}
//		if (verbosity >= 1) {
//			System.out.println();
//			System.out.println("Total iteration: " + i);
//			System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(x0), MVO.max(x0));
//			long end = System.nanoTime();
//			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
//			System.out.println();
//		}
//		if (i == maxIteration) {
//			System.out.println("Iteration Has Reached Maximum Time");
//			return -1;
//		}
//		//System.out.println("Iteration Time " + i);
//		return i;
//	}
//	
//	public int PCGNR_SSOR(double[] b, double[] x0, int verbosity) {
//		int n = b.length;
//		double[] r0 = MVO.add(b, -1.0, matU.mul(x0));
//		double[] p0 = matU.mulT(r0);
//		double[] Ap0 = new double[n];
//		int i;
//		long begin = System.nanoTime();
//		for (i = 0; i < maxIteration; i++) {
//			if (MVO.L2Norm(r0) < eps) {
//				break;
//			}
//			double temp1 = 0.0;
//			double temp2 = 0.0;
//			Arrays.fill(Ap0, 0.0);
//			matU.mulT(r0, Ap0);
//			temp1 = MVO.dot_product(Ap0, Ap0);
//			matU.mul(p0, Ap0);
//			temp2 = MVO.dot_product(Ap0, Ap0);
//			double alpha = temp1 / temp2;
//			for (int j = 0; j < n; j++) {
//				x0[j] += alpha * p0[j];
//				r0[j] -= alpha * Ap0[j];
//			}
//			
//			temp2 = temp1;
//			temp1 = 0.0;
//			//Arrays.fill(Ap0, 0.0);
//			matU.mulT(r0, Ap0);
//			temp1 = MVO.dot_product(Ap0, Ap0);
//			double beta = temp1 / temp2;
//			for (int j = 0; j < n; j++) {
//				p0[j] = Ap0[j] + beta * p0[j];
//			}
//
//			if (verbosity >= 2) {
//				System.out.println("Error = " + r0 + "    iter = " + i);
//				System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(x0), MVO.max(x0));
////				if (k == 3) {
////					System.exit(0);
////				}
//			}
////			System.out.println("Iteration " + (i + 1) + ": min\t"
////					+ MVO.min(x0) + "\tmax\t" + MVO.max(x0));
//		}
//		if (verbosity >= 1) {
//			System.out.println();
//			System.out.println("Total iteration: " + i);
//			System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(x0), MVO.max(x0));
//			long end = System.nanoTime();
//			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
//			System.out.println();
//		}
//		if (i == maxIteration) {
//			System.out.println("Iteration Has Reached Maximum Time");
//			return -1;
//		}
//		//System.out.println("Iteration Time " + i);
//		return i;
//	}
//	
//	
//	/**
//	 * �����ݶ�ƽ����(CGS)���һ�����Է�����
//	 * @param A һ�������ϡ�����
//	 * @param b �Ҷ���
//	 * @param x0 ������ʼֵ�������ص�����
//	 * @return ��������, ��������������������������-1
//	 */
//	public int CGS(double[] b, double[] x0) {
//		int n = matU.getN();
//		double[] r0 = MVO.add(b, -1.0, matU.mul(x0));
//		double[] u0 = r0.clone();
//		double[] z0 = MVO.add(r0, 1.0);
//		double[] p0 = r0.clone();
//		double[] q0 = new double[n];
//		double[] v0 = matU.mul(p0);
//		//double alpha = 0.0;
//		//double beta = 0.0;
//		int k;
//		for (k = 0; k < maxIteration; k++) {
//			if (MVO.L2Norm(r0) < eps) {
//				break;
//			}
//			
//			double temp1 = 0.0;
//			double temp2 = 0.0;
//			for(int j = 0; j < n; j++) {
//				temp1 += r0[j] * z0[j];
//				temp2 += v0[j] * z0[j];
//			}
//			double alpha = temp1 / temp2;
////			for(int i = 0; i < n; i++) {
////				q0[i] = u0[i] - alpha * v0[i];
////			}
////			double[] uq = MVO.add(u0, q0);
////			double[] Auq = A.mul(uq);
////			
////			temp2 = temp1;
////			temp1 = 0.0;
////			for(int i = 0; i < n; i++) {
////				x0[i] += alpha * uq[i];
////				r0[i] -= alpha * Auq[i];
////				temp1 += r0[i] * z0[i];
////			}
//			for(int i = 0; i < n; i++) {
//				q0[i] = u0[i] - alpha * v0[i];
//				v0[i] = u0[i] + q0[i];
//
//			}
//			
////			System.out.println("v0 = " + MVO.toString(v0));
//			//double[] uq = MVO.add(u0, q0);
//			//double[] Auq = A.mul(v0);
//			temp2 = temp1;
//			temp1 = 0.0;
//			matU.mul(v0, u0);
//			for(int i = 0; i < n; i++) {
//				x0[i] += alpha * v0[i];
//				r0[i] -= alpha * u0[i];
//				temp1 += r0[i] * z0[i];
//			}
//			double beta = temp1 / temp2;
//			for(int i = 0; i < n; i++) {
//				u0[i] = r0[i] + beta * q0[i];
//				p0[i] = u0[i] + beta *(q0[i] + beta * p0[i]);
//
//			}
//			matU.mul(p0, v0);
//			//v0 = A.mul(p0);
////			System.out.println("Iteration " + (k + 1) + ": min\t"
////					+ MVO.min(x0) + "\tmax\t" + MVO.max(x0));
//		}
//		if(k == maxIteration) {
//			System.out.println("Iteration Has Reached Maximum Time! ");
//			return -1;
//		}
//		System.out.println("Iteration Time " + k);
//		return k;
//	}
//	
//	/**
//	 * ˫�����ݶȷ�(BiCG)���һ�����Է�����
//	 * @param A һ�������ϡ�����
//	 * @param b �Ҷ���
//	 * @param x0 ������ʼֵ�������ص�����
//	 * @return ��������, ��������������������������-1
//	 */
//	public int BiCG(double[] b, double[] x0) {
//		int n = matU.getN();
//		double[] r0 = MVO.add(b, -1.0, matU.mul(x0));
//		double[] r1 = r0.clone();
//		double[] p0 = r0.clone();
//		double[] z0 = r0.clone();
//		double[] z1 = new double[n];
//		double[] pt = z0.clone();
//		int i;
//		for(i = 0; i < maxIteration; i++) {
//			if (MVO.L2Norm(r0) < eps) {
//				break;
//			}
//			double[] Ap0 = matU.mul(p0);
//			double alpha = MVO.dot_product(r0, z0)/MVO.dot_product(Ap0, pt);
//			double[] Atpt = matU.mulT(pt);
//			for(int j = 0; j < n; j++) {
//				x0[j] = x0[j] + alpha * p0[j];
//				r1[j] = r0[j] - alpha * Ap0[j];
//				z1[j] = z0[j] - alpha * Atpt[j];
//			}
//			System.out.println("z1 = " + MVO.toString(z1));
//			double beta = MVO.dot_product(r1, z1)/MVO.dot_product(r0, z0);
//			System.out.println(beta);
//			for(int j = 0; j < n; j++) {
//				p0[j] = r1[j] + beta * p0[j];
//				pt[j] = z1[j] + beta * pt[j];
//				z0[j] = z1[j];
//				r0[j] = r1[j];
//			}
//			System.out.println("Iteration " + (i + 1) + ": min\t"
//					+ MVO.min(x0) + "\tmax\t" + MVO.max(x0));
//		}
//		if (i == maxIteration) {
//			System.out.println("Iteration Has Reached Maximum Time");
//			return -1;
//		}
//		System.out.println("Iteration Time " + i);
//		return i;
//	}
//	
//	/**
//	 * ˫�����ݶ��ȶ�������������Է�����
//	 * @param A һ�������ϡ�����
//	 * @param b �Ҷ���
//	 * @param x0 ������ʼֵ�������ص�����
//	 * @return ��������, ��������������������������-1
//	 */
//	public int BiCGSTAB(double[] b, double[] x0) {
//		int n = matU.getN();
//		double[] r0 = MVO.add(b, -1.0, matU.mul(x0));
//		double[] p0 = r0.clone();
//		double[] z0 = r0.clone();
//		int i;
//		for(i = 0; i < maxIteration; i++) {
//			if(MVO.L2Norm(r0) < Constants.Er) {
//				break;
//			}
//			System.out.println("Error = " + MVO.L2Norm(r0));
//			double[] Ap0 = matU.mul(p0);
//			double temp = MVO.dot_product(r0, z0);
//			double alpha = temp/ MVO.dot_product(Ap0, z0);
//			for(int j = 0; j < n; j++) {
//				x0[j] = x0[j] + alpha * p0[j];
//				r0[j] = r0[j] - alpha * Ap0[j];
//			}
//			double[] Ar0 = matU.mul(r0);
//			double omega = MVO.dot_product(r0, Ar0)/MVO.dot_product(Ar0, Ar0);
//			for(int j = 0; j < n; j++) {
//				x0[j] = x0[j] + omega * r0[j];
//				r0[j] = r0[j] - omega * Ar0[j];
//			}
//			double beta = alpha /omega * MVO.dot_product(r0, z0) / temp;
//			for(int j = 0; j < n; j++) {
//				p0[j] = r0[j] + beta *(p0[j] - omega * Ap0[j]);
//			}
////			System.out.println("Iteration " + (i + 1) + ": min\t"
////					+ MVO.min(x0) + "\tmax\t" + MVO.max(x0));
//		}
//		
//		if (i == maxIteration) {
//			System.out.println("Iteration Has Reached Maximum Time");
//			return -1;
//		}
//		System.out.println("Iteration Time " + i);
//		return i;
//	}
	
	public void PCGSSSOR_2(Complex[] b, Complex[] x0, double omega, int verbosity) {
		int n = b.length;
		double[] bb = new double[2 * n];
		double[] xx0 = new double[2 * n];
		for(int i = 0; i < n; i++) {
			bb[i] = b[i].real;
			bb[i + n] = b[i].image;
		}
		//MVO.saveArray("AB_RHS.txt", bb);
		System.exit(1);
		NewIterSSolver solver = new NewIterSSolver(ABMatrix);
		solver.CGNR(bb, xx0, verbosity);
		//solver.CGS(bb, xx0);
		//solver.CGNE(bb, xx0, verbosity);
		//solver.PCGSSOR(bb, xx0, omega, verbosity);
		//solver.GMRES_SSOR(bb, xx0, omega, 500, verbosity);
		for(int i = 0; i < n; i++) {
			x0[i].real = xx0[i];
			x0[i].image = xx0[i + n];
		}
	}

	private CSRCMatrix matU;
	private CSRCMatrix matL;
	private Complex[] diag;
	
	private SRMatrix ABMatrix;
}
